home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / mat.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  15.3 KB  |  483 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module mat)
  13.  
  14. (COMMENT THIS IS THE MAT PACKAGE)
  15.  
  16. (DECLARE-TOP(SPECIAL PIVSIGN* *ECH* *TRI* LSOLVEFLAG $ALGEBRAIC
  17.           $MULTIPLICITIES EQUATIONS
  18.           MUL* FORMATFORM DOSIMP $DISPFLAG $RATFAC
  19.           *TB $NOLABELS ERRRJFFLAG *DET* GENVAR
  20.           XM* XN* VARLIST AX LINELABLE $LINECHAR $LINENUM SOL)
  21.      (*LEXPR $SOLVE $RAT)
  22.      (ARRAY* (NOTYPE XA* 2))
  23.      (FIXNUM TIM)
  24.      (GENPREFIX MAT))
  25.  
  26. ;;these are arrays.
  27. (DECLARE-TOP(SPECIAL *ROW* *COL* *COLINV*))
  28. ;; The array declarations of ROW, COL, and COLINV aren't having any
  29. ;; effect on the Lisp Machine.  Should be fixed somehow.
  30.  
  31. ;(DECLARE (ARRAY* (FIXNUM *ROW* 1 *COL* 1 *COLINV* 1)))
  32.  
  33. (DEFMVAR $GLOBALSOLVE NIL)
  34. (DEFMVAR $SPARSE NIL)
  35. (DEFMVAR $BACKSUBST T)
  36.  
  37. (DEFMVAR *RANK* NIL)
  38. (DEFMVAR *INV* NIL)
  39. (DEFVAR SOLVEXP NIL)
  40.  
  41. (DEFUN SOLCOEF (M *C VARL FLAG)
  42.   (PROG (CC ANSWER LEFTOVER)
  43.     (SETQ CC (CDR (RATREP* *C)))
  44.     (IF (OR (ATOM (CAR CC))
  45.         (NOT (EQUAL (CDAR CC) '(1 1)))
  46.         (NOT (EQUAL 1 (CDR CC))))
  47.         (MERROR "Unacceptable variable to SOLVE:~%~M" *C))
  48.     (SETQ ANSWER (RATREDUCE (PRODCOEF (CAR CC) (CAR M)) (CDR M)))
  49.     (IF (NOT FLAG) (RETURN ANSWER))
  50.     (SETQ LEFTOVER
  51.           (RDIS (RATPLUS M (RATTIMES (RATMINUS ANSWER) CC T))))
  52.     (IF (OR (NOT (FREEOF *C LEFTOVER))
  53.         (DEPENDSALL (RDIS ANSWER) VARL))
  54.         (ERRRJF "NON-LINEAR"))
  55.     (RETURN ANSWER)))
  56.  
  57. (DEFUN FORMX (FLAG NAM EQL VARL)
  58.   (PROG (B AX X IX J)
  59.     (SETQ VARLIST VARL)
  60.     (MAPC #'NEWVAR EQL)
  61.     (AND (NOT $ALGEBRAIC)
  62.          (ORMAPC #'ALGP VARLIST) 
  63.          (SETQ $ALGEBRAIC T))
  64.     (SET NAM (*ARRAY nil t (f1+ (SETQ XN* (LENGTH EQL)))
  65.               (f1+ (SETQ XM* (f1+ (LENGTH VARL))))))
  66.     (SETQ NAM (GET-ARRAY-POINTER NAM))
  67.     (SETQ IX 0)
  68.      LOOP1
  69.     (COND ((NULL EQL) (RETURN  VARLIST)))
  70.     (SETQ AX (CAR EQL))
  71.     (SETQ EQL (CDR EQL))
  72.     (SETQ IX (f1+ IX))
  73.     (STORE (aref NAM IX XM*) (CONST AX VARL))
  74.     (SETQ J 0)
  75.     (SETQ B VARL) (SETQ AX (CDR (RATREP* AX)))
  76.      LOOP2
  77.     (SETQ X (CAR B))
  78.     (SETQ B (CDR B))
  79.     (SETQ J (f1+ J))
  80.     (STORE (aref NAM IX J) (SOLCOEF AX X VARL FLAG))
  81.     (COND (B (GO LOOP2)))
  82.     (GO LOOP1)))
  83.  
  84. (DEFUN DEPENDSALL (EXP L)
  85.   (COND ((NULL L) NIL)
  86.     ((OR (NOT (FREEOF (CAR L) EXP)) (DEPENDSALL EXP (CDR L))) T)
  87.     (T NIL)))
  88.  
  89. (SETQ *DET* NIL *ECH* NIL *TRI* NIL)
  90.  
  91. (DEFUN PTORAT (AX M N)
  92.   (PROG (I J)
  93.     (SETQ AX (GET-ARRAY-POINTER AX))
  94.     (SETQ I (f1+ M) N (f1+ N)) 
  95.      LOOP1
  96.     (COND ((EQUAL I 1) (RETURN NIL)))
  97.     (SETQ I (f1- I) J N)
  98.      LOOP2
  99.     (COND ((EQUAL J 1) (GO LOOP1)))
  100.     (SETQ J (f1- J))
  101.     (STORE (AREF AX I J) (CONS (AREF AX I J) 1))
  102.     (GO LOOP2)))
  103.  
  104. (DEFUN MEQHK (Z)
  105.   (COND ((AND (NOT (ATOM Z)) (EQ (CAAR Z) 'MEQUAL))
  106.      (SIMPLUS (LIST '(MPLUS) (CADR Z) (LIST '(MTIMES) -1 (CADDR Z))) 1 NIL))
  107.     (T Z)))
  108.  
  109. (DEFUN CONST (E VARL)
  110.   (PROG (ZL)
  111.     (SETQ VARL (MAPCAR (FUNCTION (LAMBDA(X) (CAADR (RATREP* X)))) VARL))
  112.     (SETQ E (CDR(RATREP* E)))
  113.     (SETQ ZL (NZEROS (LENGTH VARL) NIL))
  114.     (RETURN (RATREDUCE (PCTIMES -1 (PCSUBSTY ZL VARL (CAR E)))
  115.                (PCSUBSTY ZL VARL (CDR E))))))
  116.  
  117.  
  118.  
  119. (DEFVAR *MOSESFLAG NIL)
  120.  
  121. (DEFMVAR $%RNUM 0)
  122.  
  123. (DEFMFUN MAKE-PARAM ()
  124.   (LET ((PARAM (CONCAT '$%R (SETQ $%RNUM (f1+ $%RNUM)))))
  125.     (TUCHUS $%RNUM_LIST PARAM)
  126.     PARAM))
  127.  
  128. (DEFMVAR $LINSOLVE_PARAMS T "LINSOLVE generates %Rnums")
  129.  
  130. ;(DECLARE (FIXNUM N))
  131.  
  132. (DEFUN NCDR (X N)
  133.   (NTHCDR (f1- N) X))
  134.  
  135. (DEFUN ITH (X N) (COND ((ATOM X) NIL) (T (CAR (NCDR X N)))))
  136.  
  137. ;(DECLARE (NOTYPE N))
  138.  
  139. (DEFUN POLYIZE (AX R M MUL)
  140.   (DECLARE (FIXNUM M))
  141.   (DO ((C 1 (f1+ C)) (D))
  142.       ((> C M) NIL)
  143.     (DECLARE (FIXNUM C))
  144.     (SETQ D (AREF AX R C))
  145.     (SETQ D (COND ((EQUAL MUL 1) (CAR D))
  146.           (T (PTIMES (CAR D)
  147.                  (PQUOTIENTCHK MUL (CDR D))))))
  148.     (STORE (AREF AX R C) (IF $SPARSE (CONS D 1) D))))
  149.  
  150. ; TWO-STEP FRACTION-FREE GAUSSIAN ELIMINATION ROUTINE
  151.  
  152. (DEFUN TFGELI (AX N M &AUX ($SPARSE (AND $SPARSE (OR *DET* *INV*))))
  153.   ;;$sparse is also controlling whether polyize stores polys or ratforms
  154.   (SETQ AX (GET-ARRAY-POINTER AX))
  155.   (SETQ MUL* 1)
  156.   (DO ((R 1 (f1+ R)))
  157.       ((> R N) (COND ((AND $SPARSE *DET*)(SPRDET AX N))
  158.              ((AND *INV* $SPARSE)(NEWINV AX N M))
  159.              (T (TFGELI1 AX N M))))
  160.     (DO ((C 1 (f1+ C))
  161.      (D)
  162.      (MUL 1))
  163.     ((> C M)
  164.      (AND *DET* (SETQ MUL* (PTIMES MUL* MUL)))
  165.      (POLYIZE AX R M MUL))
  166.       (COND ((EQUAL 1 (SETQ D (CDR (AREF AX R C)))) NIL)
  167.         (T (SETQ MUL (PTIMES MUL (PQUOTIENT D (PGCD MUL D)))))))))
  168.  
  169. (SETQ LSOLVEFLAG NIL)
  170.  
  171. ;; The author of the following programs is Tadatoshi Minamikawa (TM). 
  172. ;; This program is one-step fraction-free Gaussian elimination with
  173. ;; optimal pivotting.  DRB claims the hair in this program is not
  174. ;; necessary and that straightforward Gaussian elimination is sufficient,
  175. ;; for sake of future implementors.
  176.  
  177. ;; To debug, delete the comments around PRINT and BREAK statements.
  178.  
  179. (DECLARE-TOP(SPECIAL PERMSIGN A RANK DELTA NROW NVAR N M VARIABLEORDER
  180.           DEPENDENTROWS INCONSISTENTROWS L K)
  181.      ;; We could just use fortran, you know.
  182.      (FIXNUM NROW NVAR RANK I J K L M N))
  183.  
  184. (DEFUN TFGELI1 (AX N M)
  185.   (PROG (K L DELTA VARIABLEORDER INCONSISTENTROWS
  186.      DEPENDENTROWS NROW NVAR RANK PERMSIGN RESULT)
  187.     (#-cl *ARRAY #+cl cl-*array '*ROW* 'fixnum (f1+ N)) 
  188.     (#-cl *ARRAY #+cl cl-*array '*COL* 'fixnum (f1+ M))
  189.     (#-cl *ARRAY #+cl cl-*array '*COLINV* 'fixnum (f1+ M))
  190. ;    #+LISPM (FILLARRAY (FUNCTION ROW) '(0)) ;implicit in *array
  191. ;    #+LISPM (FILLARRAY (FUNCTION COL) '(0))
  192. ;    #+LISPM (FILLARRAY (FUNCTION COLINV) '(0))
  193.     (SETQ AX (GET-ARRAY-POINTER AX))
  194.         (setq *COL* (GET-ARRAY-POINTER *COL*))
  195.     (setq *ROW* (GET-ARRAY-POINTER *ROW*))
  196.     (setq *COLINV* (get-array-pointer *COLINV*))
  197.     ;; (PRINT 'ONESTEP-LIPSON-WITH-PIVOTTING)
  198.     (SETQ NROW N)
  199.     (SETQ NVAR (COND (*RANK* M) (*DET* M) (*INV* N) (*ECH* M) (*TRI* M) (T (f1- M))))
  200.     (DO ((I 1 (f1+ I))) ((> I N)) (STORE (AREF *ROW* I) I))
  201.     (DO ((I 1 (f1+ I))) ((> I M))
  202.         (STORE (AREF *COL* I) I) (STORE (AREF *COLINV* I) I))
  203.     (SETQ RESULT
  204.           (COND 
  205.         (*RANK* (FORWARD T) RANK)
  206.         (*DET* (FORWARD T)
  207.                (COND ((= NROW N) (COND (PERMSIGN  (PMINUS DELTA))
  208.                            (T DELTA)))
  209.                  (T 0)))
  210.         (*INV* (FORWARD T) (BACKWARD) (RECOVERORDER1))
  211.         (*ECH* (FORWARD NIL) (RECOVERORDER2))
  212.         (*TRI* (FORWARD NIL) (RECOVERORDER2))
  213.         (T (FORWARD T) (COND ($BACKSUBST (BACKWARD)))
  214.            (RECOVERORDER2)
  215.            (LIST DEPENDENTROWS  INCONSISTENTROWS VARIABLEORDER))))
  216.     (*REARRAY '*ROW*) (*REARRAY '*COL*) (*REARRAY '*COLINV*)
  217.     (RETURN RESULT)))
  218.  
  219. ;FORWARD ELIMINATION
  220. ;IF THE SWITCH *CPIVOT IS NIL, IT AVOIDS THE COLUMN PIVOTTING.
  221. (DEFUN FORWARD (*CPIVOT)
  222.   (SETQ DELTA 1)            ;DELTA HOLDS THE CURRENT DETERMINANT
  223.   (DO ((K 1 (f1+ K))
  224.        (NVAR NVAR)            ;PROTECTS AGAINST TEMPORARAY RESETS DONE IN PIVOT
  225.        (M M))
  226.       ((OR (> K NROW) (> K NVAR)))
  227.     (COND ((PIVOT AX K *CPIVOT) (RETURN NIL)))
  228.     ;; PIVOT IS T IF THERE IS NO MORE NON-ZERO ROW LEFT. THEN GET OUT OF THE LOOP
  229.     (DO ((I (f1+ K) (f1+ I)))
  230.     ((> I NROW))
  231.       (DO ((J (f1+ K) (f1+ J)))
  232.       ((> J M))
  233.     (STORE ( AREF AX (AREF *ROW* I) (AREF *COL* J))
  234.            (PQUOTIENT (PDIFFERENCE (PTIMES ( AREF AX (AREF *ROW* K) (AREF *COL* K))
  235.                            (AREF AX (AREF *ROW* I) (AREF *COL* J)))
  236.                        (PTIMES ( AREF AX (AREF *ROW* I) (AREF *COL* K))
  237.                            ( AREF AX (AREF *ROW* K) (AREF *COL* J))))
  238.               DELTA))))
  239.     (DO ((I (f1+ K) (f1+ I)))
  240.         ((> I NROW))
  241.       (STORE ( AREF AX (AREF *ROW* I) (AREF *COL* K))
  242.          0))
  243.     (SETQ DELTA ( AREF AX (AREF *ROW* K) (AREF *COL* K))))
  244.   ;; UNDOES COLUMN HACK IN PIVOT.
  245.   (OR *CPIVOT (DO ((I 1 (f1+ I))) ((> I M)) (STORE (AREF *COL* I) I)))
  246.   (SETQ RANK (MIN NROW NVAR)))
  247.  
  248. ; BACKWARD SUBSTITUTION
  249. (DEFUN BACKWARD ()
  250.   (DO ((I (f1- RANK) (f1- I)))
  251.       ((< I 1))
  252.     (DO ((L (f1+ RANK) (f1+ L)))
  253.     ((> L M))
  254.       (STORE (AREF AX (AREF *ROW* I) (AREF *COL* L))
  255.          (PQUOTIENT (PDIFFERENCE
  256.               (PTIMES (AREF AX (AREF *ROW* I) (AREF *COL* L))
  257.                   (AREF AX (AREF *ROW* RANK) (AREF *COL* RANK)))
  258.               (DO ((J (f1+ I) (f1+ J)) (SUM 0))
  259.                   ((> J RANK) SUM)
  260.                 (SETQ SUM (PPLUS SUM (PTIMES (AREF AX (AREF *ROW* I) (AREF *COL* J))
  261.                              (AREF AX (AREF *ROW* J) (AREF *COL* L)))))))
  262.             (AREF AX (AREF *ROW* I) (AREF *COL* I)))))
  263.     (DO ((L (f1+ I) (f1+ L)))
  264.     ((> L RANK))
  265.       (STORE (AREF AX (AREF *ROW* I) (AREF *COL* L)) 0)))
  266.   ;; PUT DELTA INTO THE DIAGONAL MATRIX
  267.   (SETQ DELTA (AREF AX (AREF *ROW* RANK) (AREF *COL* RANK)))
  268.   (DO ((I 1 (f1+ I)))
  269.       ((> I RANK))
  270.     (STORE (AREF AX (AREF *ROW* I) (AREF *COL* I)) DELTA)))
  271.  
  272. ;RECOVER THE ORDER OF ROWS AND COLUMNS.
  273.  
  274. (DEFUN RECOVERORDER1 ()
  275.   ;;(PRINT 'REARRANGE)
  276.   (DO ((I NVAR (f1- I)))
  277.       ((= I 0))
  278.     (SETQ VARIABLEORDER (CONS I VARIABLEORDER)))
  279.   (DO ((I (f1+ RANK) (f1+ I)))
  280.       ((> I N))
  281.     (COND ((EQUAL (AREF AX (AREF *ROW* I) (AREF *COL* M)) 0) 
  282.        (SETQ DEPENDENTROWS (CONS (AREF *ROW* I) DEPENDENTROWS)))
  283.       (T (SETQ INCONSISTENTROWS (CONS (AREF *ROW* I) INCONSISTENTROWS)))))
  284.   (DO ((I 1 (f1+ I)))
  285.       ((> I N))
  286.     (COND ((NOT (= (AREF *ROW* (AREF *COLINV* I)) I))
  287.        (PROG ()
  288.          (MOVEROW AX N M I 0)
  289.          (SETQ L I)
  290.           LOOP
  291.          (SETQ K (AREF *ROW* (AREF *COLINV* L)))
  292.          (STORE (AREF *ROW* (AREF *COLINV* L)) L)
  293.          (COND ((= K I) (MOVEROW AX N M 0 L))
  294.                (T (MOVEROW AX N M K L)
  295.               (SETQ L K)
  296.               (GO LOOP))))))))
  297.  
  298. (DEFUN RECOVERORDER2 ()
  299.   (DO ((I NVAR (f1- I)))
  300.       ((= I 0))
  301.     (SETQ VARIABLEORDER (CONS (AREF *COL* I) VARIABLEORDER)))
  302.   (DO ((I (f1+ RANK) (f1+ I)))
  303.       ((> I N))
  304.     (COND ((EQUAL (AREF AX (AREF *ROW* I) (AREF *COL* M)) 0)
  305.        (SETQ DEPENDENTROWS (CONS (AREF *ROW* I) DEPENDENTROWS)))
  306.       (T (SETQ INCONSISTENTROWS (CONS (AREF *ROW* I) INCONSISTENTROWS)))))
  307.   (DO ((I 1 (f1+ I)))
  308.       ((> I N))
  309.     (COND ((NOT (= (AREF *ROW* I) I))
  310.        (PROG ()
  311.          (MOVEROW AX N M I 0)
  312.          (SETQ L I)
  313.           LOOP
  314.          (SETQ K (AREF *ROW* L))
  315.          (STORE (AREF *ROW* L) L)
  316.          (COND ((= K I) (MOVEROW AX N M 0 L))
  317.                (T (MOVEROW AX N M K L)
  318.               (SETQ L K)
  319.               (GO LOOP)))))))
  320.   (DO ((I 1 (f1+ I)))
  321.       ((> I NVAR))
  322.     (COND ((NOT (= (AREF *COL* I) I))
  323.        (PROG ()
  324.          (MOVECOL AX N M I 0)
  325.          (SETQ L I)
  326.           LOOP2
  327.          (SETQ K (AREF *COL* L))
  328.          (STORE (AREF *COL* L) L)
  329.          (COND ((= K I) (MOVECOL AX N M 0 L))
  330.                (T (MOVECOL AX N M K L)
  331.               (SETQ L K)
  332.               (GO LOOP2))))))))
  333.  
  334. ;THIS PROGRAM IS USED IN REARRANGEMENT
  335. (DEFUN MOVEROW (AX N M I J)
  336.   (DO ((K 1 (f1+ K))) ((> K M))
  337.     (STORE (AREF AX J K) (AREF AX I K))))
  338.  
  339. (DEFUN MOVECOL (AX N M I J)
  340.   (DO ((K 1 (f1+ K))) ((> K N))
  341.     (STORE (AREF AX K J) (AREF AX K I))))
  342.  
  343. ;COMPLEXITY IS DEFINED AS FOLLOWS
  344. ; COMPLEXITY(0)=0
  345. ; COMPLEXITY(CONSTANT)=1
  346. ; COMPLEXITY(POLYNOMIAL)=1+SUM(COMPLEXITY(C(N))+COMPLEXITY(E(N)), FOR N=0,1 ...M)
  347. ; WHERE POLYNOMIAL IS OF THE FORM
  348. ;    SUM(C(N)*X^E(N), FOR N=0,1 ... M)     X IS THE VARIABLE
  349.  
  350. (DEFUN COMPLEXITY (EXP)
  351.   (COND ((NULL EXP) 0)
  352.     ((EQUAL EXP 0) 0)
  353.     ((ATOM  EXP) 1)
  354.     (T (PLUS (COMPLEXITY (CAR EXP)) (COMPLEXITY (CDR EXP))))))
  355.  
  356. (DEFUN COMPLEXITY/ROW (AX I J1 J2)
  357.   (DO ((J J1 (f1+ J)) (SUM 0))
  358.       ((> J J2) SUM)
  359.     (SETQ SUM (PLUS SUM (COMPLEXITY (AREF AX (AREF *ROW* I) (AREF *COL* J)))))))
  360.  
  361. (DEFUN COMPLEXITY/COL (AX J I1 I2)
  362.   (DO ((I I1 (f1+ I)) (SUM 0))
  363.       ((> I I2) SUM)
  364.     (SETQ SUM (PLUS SUM (COMPLEXITY (AREF AX (AREF *ROW* I) (AREF *COL* J)))))))
  365.  
  366. (DEFUN ZEROP/ROW (AX I J1 J2)
  367.    (DO ((J J1 (f1+ J)))
  368.        ((> J J2) T)
  369.       (COND ((NOT (EQUAL (AREF AX (AREF *ROW* I) (AREF *COL* J)) 0)) (RETURN NIL)))))
  370.  
  371. ;PIVOTTING ALGORITHM
  372. (DEFUN PIVOT (AX K *CPIVOT)
  373.   (PROG ( ROW/OPTIMAL COL/OPTIMAL COMPLEXITY/I/MIN COMPLEXITY/J/MIN
  374.      COMPLEXITY/I COMPLEXITY/J COMPLEXITY/DET COMPLEXITY/DET/MIN DUMMY)
  375.     (SETQ ROW/OPTIMAL K COMPLEXITY/I/MIN 1000000. COMPLEXITY/J/MIN 1000000.)
  376.     ;;TEST THE SINGULARITY
  377.     (COND ((DO ((I K (f1+ I)) (ISALLZERO T))
  378.            ((> I NROW) ISALLZERO)
  379.          LOOP (COND ((ZEROP/ROW AX I K NVAR)
  380.                  (COND (*INV* (MERROR "Singular"))
  381.                    (T (EXCHANGEROW I NROW)
  382.                       (SETQ NROW (f1- NROW))))
  383.                  (COND ((NOT (> I NROW)) (GO LOOP))))
  384.                 (T (SETQ ISALLZERO NIL))))
  385.            (RETURN T)))
  386.  
  387.     ;; FIND AN OPTIMAL ROW
  388.     ;; IF *CPIVOT IS NIL, (AX I K) WHICH IS TO BE THE PIVOT MUST BE NONZERO.
  389.     ;; BUT IF *CPIVOT IS T, IT IS UNNECESSARY BECAUSE WE CAN DO THE COLUMN PIVOT.
  390.      FINDROW
  391.     (DO ((I K (f1+ I)))
  392.         ((> I NROW))
  393.       (COND ((OR *CPIVOT (NOT (EQUAL (AREF AX (AREF *ROW* I) (AREF *COL* K)) 0)))
  394.          (COND ((> COMPLEXITY/I/MIN
  395.                (SETQ COMPLEXITY/I (COMPLEXITY/ROW AX I K M)))
  396.             (SETQ ROW/OPTIMAL I COMPLEXITY/I/MIN COMPLEXITY/I))))))
  397.     ;; EXCHANGE THE ROWS K AND ROW/OPTIMAL
  398.     (EXCHANGEROW K ROW/OPTIMAL)
  399.  
  400.     ;; IF THE FLAG *CPIVOT IS NIL, THE FOLLOWING STEPS ARE NOT EXECUTED.
  401.     ;; THIS TREATMENT WAS DONE FOR THE LSA AND ECHELONTHINGS WHICH ARE NOT
  402.     ;; HAPPY WITH THE COLUMN OPERATIONS.
  403.     (COND ((NULL *CPIVOT)
  404.            (COND ((NOT (EQUAL (AREF AX (AREF *ROW* K) (AREF *COL* K)) 0))
  405.               (RETURN NIL))
  406.              (T (DO ((I K (f1+ I))) ((= I NVAR))
  407.                 (STORE (AREF *COL* I) (AREF *COL* (f1+ I))))
  408.             (SETQ NVAR (f1- NVAR) M (f1- M))
  409.             (GO FINDROW)))))
  410.  
  411.     ;;STEP3 ... FIND THE OPTIMAL COLUMN
  412.     (SETQ COL/OPTIMAL 0
  413.           COMPLEXITY/DET/MIN 1000000.
  414.           COMPLEXITY/J/MIN 1000000.)
  415.  
  416.     (DO ((J K (f1+ J)))
  417.         ((> J NVAR))
  418.       (COND ((NOT (EQUAL (AREF AX (AREF *ROW* K) (AREF *COL* J)) 0))
  419.          (COND ((> COMPLEXITY/DET/MIN
  420.                (SETQ COMPLEXITY/DET
  421.                  (COMPLEXITY (AREF AX (AREF *ROW* K) (AREF *COL* J)))))
  422.             (SETQ COL/OPTIMAL J
  423.                   COMPLEXITY/DET/MIN COMPLEXITY/DET
  424.                   COMPLEXITY/J/MIN (COMPLEXITY/COL AX J (f1+ K) N)))
  425.                ((EQUAL COMPLEXITY/DET/MIN COMPLEXITY/DET)
  426.             (COND ((> COMPLEXITY/J/MIN
  427.                   (SETQ COMPLEXITY/J
  428.                     (COMPLEXITY/COL AX J (f1+ K) N)))
  429.                    (SETQ COL/OPTIMAL J
  430.                      COMPLEXITY/DET/MIN COMPLEXITY/DET
  431.                      COMPLEXITY/J/MIN COMPLEXITY/J))))))))
  432.  
  433.     ;(COND ((ZEROP COL/OPTIMAL) (COMMENT (PRINT '"SINGULAR!"))))
  434.  
  435.     ;; EXCHANGE THE COLUMNS K AND COL/OPTIMAL
  436.     (EXCHANGECOL  K COL/OPTIMAL)
  437.     (SETQ DUMMY (AREF *COLINV* (AREF *COL* K)))
  438.     (STORE (AREF *COLINV* (AREF *COL* K)) (AREF *COLINV* (AREF *COL* COL/OPTIMAL)))
  439.     (STORE (AREF *COLINV* (AREF *COL* COL/OPTIMAL)) DUMMY)
  440.     (RETURN NIL)))
  441.  
  442. (DEFUN EXCHANGEROW (I J)
  443.   (PROG (DUMMY)
  444.     (SETQ DUMMY (AREF *ROW* I))
  445.     (STORE (AREF *ROW* I) (AREF *ROW* J))
  446.     (STORE (AREF *ROW* J) DUMMY) 
  447.     (COND ((= I J) (RETURN NIL))
  448.           (T (SETQ PERMSIGN (NOT PERMSIGN))))))
  449.  
  450. (DEFUN EXCHANGECOL (I J)
  451.   (PROG (DUMMY)
  452.     (SETQ DUMMY (AREF *COL* I))
  453.     (STORE (AREF *COL* I) (AREF *COL* J))
  454.     (STORE (AREF *COL* J) DUMMY)
  455.     (COND ((= I J) (RETURN NIL))
  456.           (T (SETQ PERMSIGN (NOT PERMSIGN))))))
  457.  
  458. ;; Displays list of solutions.
  459.  
  460. (DEFUN SOLVE2 (Llist)
  461.   (SETQ $MULTIPLICITIES NIL)
  462.   (MAP2C #'(LAMBDA (EQUATN MULTIPL) 
  463.          (SETQ EQUATIONS
  464.            (NCONC EQUATIONS (LIST (DISPLINE EQUATN))))
  465.          (PUSH MULTIPL $MULTIPLICITIES)
  466.          (IF (AND (> MULTIPL 1) $DISPFLAG)
  467.          (MTELL "Multiplicity ~A~%" MULTIPL)))
  468.      Llist)
  469.   (SETQ $MULTIPLICITIES (CONS '(MLIST SIMP) (NREVERSE $MULTIPLICITIES))))
  470.  
  471. ;; Displays an expression and returns its linelabel.
  472.  
  473. (DEFMFUN DISPLINE (EXP)
  474.  (LET ($NOLABELS (TIM 0)) 
  475.       (ELABEL EXP)
  476.       (COND ($DISPFLAG (REMPROP LINELABLE 'NODISP) 
  477.                (SETQ TIM (RUNTIME))
  478.                (MTERPRI)
  479.                (DISPLA (LIST '(MLABLE) LINELABLE EXP))
  480.                (TIMEORG TIM))
  481.         (T (PUTPROP LINELABLE T 'NODISP)))
  482.       LINELABLE))
  483.